Optimal reduction and conversion of range-difference measurements for positioning

For positioning an object with m references, there are m−1 linearly independent range differences and measuring them is essential. However, none of m(m−1) possible range differences should be considered redundant unless their measurements are free of noise and locations of the references are exactly known. From all available range-difference measurements, m range measurements are obtained for positioning based on the least squares principle. The problem formulation treats missing and weighted range-difference measurements simultaneously. The exact relationships among several formulations of least squares positioning are established. A numerical example illustrates the results.


Introduction
Positioning an energy-emitting or reflecting object is an intensively studied topic due to its importance in wide applications [1]. Often it is based on the principle of time differences of arrival (TDOA), which means use of indirect measurements of range differences between the object at an unknown location and references at known locations. Positioning with TDOA measurements is different but closely related to that using time of arrival (TOA) measurements with or without a bias. With all possible range measurements between each object pair, positioning multiple objects is possible [2].
Positioning an nD object with m references requires m > n and up to m(m−1) range differences can be formed, but only m−1 of them are linearly independent. Due to noise effects, measurements of all available range differences should be used for positioning in applications. This work shows how to combine all available TODA measurements to form m TOA measurements for least squares positioning. Examinations of least squares criteria of several types of TOA and TDOA measurement equations establish equivalent and other exact relationships among these positioning formulations. The cases of positioning with missing and weighted TDOA measurements are treated inclusively in this work.
Investigation on the underlying problem is important from both theoretical and practical viewpoints. This kind of study answers the question of whether or not different TODA and TOA formulations are equivalent for positioning, and provides simplified equations for algorithm development and implementation in applications. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Related work
For positioning an nD object with m known references, most methods have used m TOA equations or m−1 TODA equations for positioning, and normally the minimum number of references m = n + 1 is assumed. In the majority of TDOA methods, measurements of the remaining possible (m−1) 2 range differences are unused or assumed to be unavailable.
In an early study [3], a large number of TDOA measurements with the minimum number of references were combined to form TODA triads for improvement on positioning. Optimality of the combination was not addressed nevertheless.
The problem of TDOA denoising [4] is to find a range-measurement vector for generation of an ideally structured TDOA measurement matrix closest to the original noise-corrupted same matrix by least squares. This problem is related but not equivalent to the TDOA positioning problem directly addressed in the current study.
The problem formulation in the current study avoids the assumptions on skew symmetry of the noise-corrupted TDOA measurement matrix, and on Gaussian distribution of the noise [4,5]. This allows a more general coverage of noise conditions and consideration of up to m(m −1) TDOA measurements rather than half of them. Normally missing and weighted TDOA measurements are treated separately, for instance, in [4,6], but simultaneously in the current study.
The focus of this work is on optimal conversion of range-difference equations rather than solving them. This is because, assuming exact solvability, closed-form solutions are known for positioning with m biased TOA measurements [7][8][9][10][11][12], and with m − 1 TDOA measurements [13][14][15]. As well known, m TOA measurements can be trivially converted to m − 1 TDOA measurements, although optimality of such a conversion is unclear. In real applications, closed-form solutions offer fine approximations, and can also be used to initiate an iterative algorithm for improving the solutions. These methods can be applied to the m TOA or further m − 1 TDOA measurement equations converted from possible m(m − 1) range-difference equations studied in the current work.

Notations
All considered quantities are real numbers. Scalars are lowercase letters, (column) vectors and matrices are boldfaced lowercase and uppercase letters respectively. Set {a i } contains elements with a known number, and they can form a vector a = [a i ]. The vector of ones is denoted by e. A = [a i,j ] is a matrix of a known size with a ij being its element in the ith row and jth column. Diagonal matrix D a has elements of a on its diagonal, and D e is the identity matrix I. A 0 , trA, rankA and A + are the transpose, trace, rank and Moore-Penrose inverse of A respectively. The norm of a is jaj ¼ ffi ffi ffi ffi ffi ffi a 0 a p , and that of A is jAj ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi trðA 0 AÞ p . Denoted by A � B is the entry-wise multiplication, namely Hadamard product of the two matrices of the same size. A ¼ A 1 2 ð Þ 0 A 1 2 is a positive semi-definite matrix decomposed by its square root (matrix). A = USV 0 is the singular value decomposition of A with U 0 U = I, V 0 V = I, and diagonal matrix S consisting of (nonnegative) singular values of A. Denoted by arg min x f(x) is the argument of the minimum of a scalar function, namely x minimizing f(x). Denoted by n i � N ð� n; s 2 n Þ is a random variable n i satisfying the Gaussian distribution with mean � n and variance s 2 n . Similarly, n � N ð� n; S n Þ stands for a Gaussian random vector with mean � n and variance matrix S n .

Problem formulation
Denote the matrices of range differences and their measurements by R = [r ij ] and T = [τ ij ], respectively. On the TDOA principle, noisy measurements {τ ij } of range differences {r ij } are described by scalar equations where n ij is a random variable with zero mean, m the number of references, and r ij = r i − r j the difference of ranges r k = |p − p k | for k = i, j, from an unknown object p to known references p i and p j . If p is of dimension n, m > n is required. To have a unique p in the noise-free case, {p i } are assumed to be non-coplanar, namely they are not located in an (n − 1)D linear subspace.
Likely not all m(m − 1) measurements in (1) are available in applications even under the assumption that, without loss of generality, all m references have been used in generation of the measurements. Available measurements may also be weighted according to a priori knowledge of noise statistics. To consider cases of missing and weighted measurements simultaneously, define a masking matrix as where weight w ij > 0, and in the case of non-weighting w ij = 1. Weights {w ij } in (2) could be chosen as the components of the inverse variance matrix of noise {n ij } in (1). This resembles the treatment of measurement noise in the Kalman filtering. However, the problem considered in this study is not the problem of tracking a moving target because, if any, dynamics of object p is not considered in the current study. Hence, positioning an object based on the equivalent range equations is generally not the minimum variance estimation intended with a Kalman filter.
Based on all available range-difference measurements and possibly also with weighting, positioning an object is to find a least squares solution of p to the matrix equation The objective of this study is to convert (3) which may have up to m(m − 1) range-difference equations to m range equations. The conversion is optimal in the sense of least squares.
In the case where no range-difference measurement is weighted or missing, the scalar equations in (1) are identical to the matrix equation in (3) except that the noise terms in the former are set to be zero in the latter. In general, (3) is a compact notation of (1) by setting the unknown noise to be zero but with simultaneous consideration of weighted and missing measurements. Clearly, (3) does not normally have an exact solution for p, and hence an estimation of p is sought in respect to least squares of (3).

Linear dependence of {r ij } and properties of E
Range differences {r ij } are clearly related to each other, and linear independence of a subset of them is defined conventionally.
Definition 1 fr i l j l g for i l , j l 2 {1, 2, . . ., m} and l = 1, 2, . . ., k with an arbitrary integer k > 0, are said to be linearly independent from each other if P k l¼1 a l r i l j l ¼ 0 implies coefficient α l = 0 for all l.
Masking matrix E is non-negative, namely none of its components is negative. Also, E 6 ¼ 0 because at least two range-difference measurements are available under the necessary condition m > n for unique positioning of an nD object with m references. If every reference has been used in the generation of TDOA measurements, for all i, the ith row and column of E cannot be simultaneously zero. This amounts to, for at least one j, The case e ij + e ji = 0 for a particular i and all j, corresponds to non-use of the ith reference, which can be handled by dropping p i and reducing number m by one in (1). As implied in (4), measurements of m − 1 linearly independent range differences are automatically available in (3).

Pseudo range-measurement vector τ and properties of companion matrix � E
Define a companion matrix of E as and a pseudo range measurement vector as where � E þ is the Moore-Penrose inverse of � E, and e the vector of ones. In the special case where all m(m − 1) measurements in (1) are available and no weighting is applied to them, it is ready to obtain the simplifications E = ee 0 − I, E is obviously symmetric, and � E 6 ¼ 0 due to E 6 ¼ 0. To explore its properties, some basic definitions related to matrix irreducibility are needed. These properties are important for reduction and conversion of the weighted range difference matrix equation in (3).

Equivalence and optimality of range and range-difference equations
In terms of an arbitrary vector x of dimension m, define a matrix as which has the same structure as R, in fact R = T r with range r = [r i ]. For least squares positioning, exact relationships among (3) and the following three equations with � r ¼ e 0 ðτ À rÞ=m, are shown in the next theorem. Theorem 2 There are two least squares positioning equivalences: arg min For arbitrary p, r, and E, the following relations hold The significance of Theorem 2 lies in establishment of the equivalence of matrix Eq (3) and that in (8) to the vector equations in (8) respectively through (9) and (10) for least squares positioning. Although further equivalence between (9) and (10) cannot be established, (11) implies that if a p diminishes jτ À r À � rej considerably, it is a superb approximation of (9), and confirms the supremacy of τ ¼ r þ � re over τ ¼ r þ re and τ ¼ r for determination of p by the least squares principle.

Corollary 1
The denoising problem has the general solution with an arbitrary r. The result presented in Corollary 1 was first obtained in [4], and now given without imposing any particular assumptions on T in (1) and E in (4). Theorem 2 and Corollary 1 indicate exactly the relationship between the positioning and denoising problems. Basically, for positioning p ¼ argmin p jE � ðT À T r Þj 2 , while for denoising, x ¼ argmin x jE � ðT À T x Þj 2 . In general, jE � ðT À T r Þj 2 6 ¼ jE � ðT À T τÀ re Þj 2 , and the equality holds if p satisfies τ ¼ r þ re for some r. Note that normally τ ¼ r þ re is not exactly solvable for p and r.
It is well known that biased TOA equations are often described by τ ¼ r þ re with r representing the clock bias between the transmitter and receiver. Interestingly, in τ ¼ r þ � re, � r is specified as the average of the difference between {τ i } and {r i }. As implied in the proof of Theorem 2, � r is actually the least squares solution of r to τ ¼ r þ re.

Theoretical verification
Some primary results on matrix irreducibility are needed for proving Theorem 1.

Proof of Theorem 1
Direct calculations give where � E e ¼ ½e 2 ij �, and D � e is the diagonal matrix formed by � e ¼ ½� e i � with Obviously, � E is symmetric, and diagonally dominant, which is a). It is also at least positive semi-definite which is part of b) due to following from some simple properties of Hadamard products [16], such as and identities for arbitrary vectors x and y, and arbitrary matrices A, B and C, all with compatible dimensions. This verifies c) due to the existence of decomposition � It also leads to b) because of rank deficiency of � E in view of � Ee ¼ 0 from (13). Moreover, � τ 0 e ¼ 0 follows from the skew symmetry of E � E � T À ðE � E � TÞ 0 , which implies rank½ � E; � τ � < m. A deductive verification of irreducibility of � E and rank � E ¼ m À 1 in the following completes the proof of d) and e).
Suppose � E k is irreducible and rank � E k ¼ k À 1 for k > 2. Trivially, irreducibility of � E k implies the same of � E k þ D � e kþ1 . In view of (14), � e kþ1 6 ¼ 0 and � e kþ1 6 ¼ 0, which ensures irreducibility of � E kþ1 according to Lemma 2. Moreover, � E k þ D � e kþ1 is irreducibly diagonally dominant, and hence non-singular according to Lemma 1. This verifies rank � E kþ1 � k, where the inequality cannot hold nevertheless because � E kþ1 e ¼ 0 follows from (13).

Proof of Theorem 2
By definition, the left side of (9) specifies least squares solutions to (3). Using (6), (16), (17), D e = I, and some simple properties of the trace of matrix products, the following is obtained Since the first two terms in (23) are independent of p, the equivalence in (9) is then proved.
À r À reÞ for arbitrary p and r. By setting r ¼ � r and noticing Ab| � |A||b| for arbitrary A and b, the first inequality in (11) is verified. Setting partial differentiation @jτ À r À rej 2 =@r ¼ 0 leads to r ¼ � r, which implies jτ À r À rej � jτ À r À � rej for arbitrary p and r, and hence confirms the second inequality. From this, the third inequality follows immediately.

Proof of Corollary 1
Noting R = T r and from (23) with T x replacing R, it is ready to have arg min The general least squares solution of x to � E 1 2 ðτ À xÞ ¼ 0 is where y is arbitrary but subject to � E 1 2 y ¼ 0. Considering singular value decomposition � Recalling � E 1 2 e ¼ 0 and rank � E 1 2 ¼ m À 1, y must be parallel to E, which leads to the general expression y ¼ À re with r being an arbitrary scalar.

Illustrative example
A numerical example of 3D positioning is used to illustrate the developed results. When used for illustrating effects of different sets of range and range-difference equations on positioning, simulated datasets are considered most effective than practical datasets. This is because noise levels for range-difference measurements and inaccuracy of the reference locations could be easily set and examined in numerical examples. It is however not the case in a real setup where inaccuracies of the measurements and reference locations are coupled with the uncertainty of the object location.
Let four references be located precisely at ½ � p 1 ; � p 2 ; � p 3 ; � p 4 � ¼ ½ 0; I �, but inexactly known as Measurements of range differences are produced according to τ ij = r ij + n ij in (1) with r ij = r i −r j and r k ¼ jp À � p k j for k = i, j and noise variable n ij � N ð0; s 2 n Þ. An object is with its coordinates randomly generated within range (1, 10) as p ¼ 9:8115 6:0293 1:1923 � 0 : ð25Þ � The standard deviations are set as σ p = 0.1% and σ n = max{|r ij |} × 1%. For instance, one simulation run produced reference and measurement matrices as : As expected, since the noise levels are low, deviations of {p i } from f� p i g are insignificant, and also due to r ij = −r ji , T is approximately skew symmetric. To find the object position, the unconstrained multivariable minimization algorithm fminsearch in MATLAB has been used. In numerical minimization, referring to Theorem 2, the squared forms of the criteria aÞ jE�ðTÀ RÞj; bÞ j � E 1 2ðτ À rÞj; cÞ jT τ À Rj; dÞ jτÀ rÀ � rej; eÞ jτÀ rÀ rej fÞ jτÀ rj were used, and the initial estimate of (p, r) was taken as ðp=2; � r=2Þ in each simulation run. Fig 1 indicates the geometric setup for positioning a 3D object with four references. Table 1 shows the estimates of the object position in 50 runs of simulations under the noise conditions stated above for the reference locations and rage-difference measurements. The estimates using the squared criteria (26a) and (26b) are very close to each other, while those using criteria (26b) to (26e) are indistinguishable from each other. Due to r = − 11.81 on average over the 50 simulation runs, the estimates using the squared criteria (26f) are too poor to be useful. This value of r should not be interpreted as a clock bias because the generation of the ranges and their noisy measurements had not introduced an offset in each run of the simulations. In fact, even if introduced, a bias in TOA measurements cannot be recovered from TDOA measurements.
The best and worst estimates are determined with respect to jp Àpj which cannot be evaluated in real applications nevertheless. Worst cases of the randomly generated reference locations and range difference measurements ought to be responsible for worst estimates of p in Table 1. This is because the evaluations of criteria (26b) to (26e) have generated insignificant values at level 10 −7 which corresponds to level 10 −14 produced by least squares. As expected, if no noise is added to the reference locations and range-difference measurements, all estimates, In Table 2, three averaged values of each minimization criterion and estimation error correspond to the above three cases. On average, use of more measurements is shown to have better estimations. This indicates, as expected, that all available range-difference measurements should be used for positioning. Use of different squared criteria in (26), except for (26f), in minimization has produced estimations of the object position with similar or identical

Concluding remarks
Given m references, following the procedure in [18] on the basis of processing individually received signals, up to m(m − 1)/2 TDOA measurements could be made available. The procedure in [19] on the basis of processing each pair of received signals could produce up to m(m − 1) TDOA measurements. The current work has proposed a general method for use of these multiple TDOA measurements for positioning. While weighted least squares positioning has been considered in this work, it has not addressed issues of selection of weighting coefficients. As widely used, for instance in [5,6] among others, an obvious choice of weightings is the inverse variances of the measurement noise components, which could be obtained from a processor of generating TDOA measurements such as those in [18,19].
Minimizing the difference between the noise-corrupted TDOA measurement matrix and the well-structured matrix formed by a range-measurement vector is the denoising problem explored in [4]. The current work directly addresses the problem of minimizing TDOA equations with respect to the object location and automatically obtains the range red-measurement vector. The current focus is on positioning in a general setting with a simultaneous treatment of missing and weighted TDOA measurements. The current work has however not considered the issue of eliminating outlier measurements as examined in [4] and comprehensively explored in [20].
A numerical example has been used to illustrate the theoretical results presented in this paper. To evaluate effectiveness of the proposed method in real applications, an experiment could be designed. It would need a high-precision positioning system for referencing, and apply the method to a low-precision TDOA dataset. A practical application of the results in this paper could be in wireless sensor networks [21] with massive low-cost miniature sensors often randomly deployed in a geographical area, where a sensor could be localized by using a huge number of TDOA measurements from references including already localized sensor nodes [22].